Long- Range Order and Interactions of Macroscopic Objects in Polar Liquids. 
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We develop a phenomenological vector model of polar liquids capable to describe aqueous inter- 
actions of macroscopic bodies. It is shown that a strong, long-range and orientationally dependent 
interaction between macroscopic objects appears as a result of competition between short-range 
(hydrogen bonding) and the long-range dipole-dipole interactions of the solvent molecules. Spon- 
taneous polarization of molecular dipoles next to a hydrophobic boundaries leads to formation of 
globally ordered network of hydrogen-bonded molecules with ferroelectric properties. The proposed 
vector model naturally describes topological excitations on the solute boundaries and can be used 
to explain the hydrogen bonds networks and order-disorder phase transitions in the hydration water 
layer. 



INTRODUCTION. 



Interactions of macroscopic bodies in aqueous environ- 
ments is a fundamentally important problem in physics, 
chemistry, structural biology and in silico drug design. 
The main theoretical difficulty arises both from strong 
long-range dipole-dipole interactions of molecular dipole 
moments and complicated nature of short-range interac- 
tion between water molecules. For practical purposes, 
e.g. calculation of binding free energy (inhibition con- 
stant) of a protein-ligand complex, protonation states 
predictions for biomacromolecules, physics of membranes 
etc, the ultimate way to quantitatively account for sol- 
vation effects is the Free Energy Perturbation method 
(FEP) based on a direct modelling of both the solute and 
the solvent molecules, molecular dynamics (MD). One of 
the main advantages of MD is a possibility to include the 
surrounding water molecules into the calculation directly 
(simulations with explicit water) . In fact this is so far one 
of the most accurate ways to compute the solvation and 
the binding energies Q. The limitation of the method 
comes together with its strength: to do it one should 
compute coordinates and velocities of very large num- 
ber of atoms with time step ts ~ 10 _15 s. To allow the 
aqueous environment to relax fully, the averaging needs 
to be performed on a sufficiently long time span to in- 
clude at least the life-time of hydrogen bonds (~ 10 _9 s), 
and even the much larger relaxation and rearranging time 
of macroscopic water clusters (up to ~ lCP 6 s, see discus- 
sion below) . For modern computers it is possible to follow 
the evolution of macromolecules in water environment up 
to times ~ 10 _7 s, which may not be sufficient to relax 
the surrounding water molecules in some situations. Al- 
though such calculations can be very accurate, realistic 
applications require enormous computational resources. 



Polar liquids such as water are characterized by large 
values of static dielectric constants e 3> 1 . That is why 
to a large extent solvation theory can be reduced to a 
macroscopic electrostatics with a solvent being modelled 
as a high-e media and the solute treated as being emerged 
in a low-e cavity. This approach dates back to early 



Born's papers and can be very successful for quantita- 
tive predictions of solvation energies of small molecules 
or biomolecules interactions 0, HI- On other hand the 
interactions of sufficiently small neutral objects is domi- 
nated by hydrophobic effect, which according to 0. 
originates from short-range interactions of the solvent 
molecules. Density functional models 0. 0- @ success- 
fully explain the hydrophobic interactions on molecular 
scales, though inclusion of electrostatic interactions ap- 
pears to be a difficult task. Moreover, understanding of 
a number of important observations revealed from MD 
simulations must require complete inclusion of long-range 
interactions and thus a more sophisticated approach. 
Among the examples are the arrangement of the water 
molecules dipole moments parallel to the hydrophobic 
surface the vortex- like structures of molecular 

dipole moments networks and dipole-bridges near and 
between the solvated objects [Tri ]. 

In this Letter we introduce a continuous vector model 
of a polar liquid capable both the short- and the long- 
range features of a polar liquid in a single theoreti- 
cal framework. The model can be applied to aque- 
ous interactions of macroscopic bodies of various shapes 
and charges. It is shown that the competition between 
the short range (hydrogen bonding) and the long-range 
dipole-dipole interactions of the solvent molecules leads 
to appearance of strong, long range and orientationally 
dependent interactions between macroscopic objects. A 
polar liquid itself is characterized by a complicated fluc- 
tuating thermal state, ordered at sufficiently short scales 
within a single domain, and completely disordered at 
larger distances. This physical picture has far reaching 
consequences, especially at solvent-solute surfaces, where 
the ferroelectric film of solvent molecules is formed. With 
temperature increasing fluctuations of molecular dipole 
moments can lead to topological phase transitions de- 
pending on hydrophobic properties of the interface. We 
argue that the dynamics of macroscopic topological ex- 
citations on a surface of the liquid may be a physical 
picture behind the percolation transition of the hydrogen 
bonding network observed in MD simulations [l^lLall3l |. 

The manuscript is organized as follows. After the in- 
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traduction in Section 2 we discuss the basic assumptions 
and formulate the vector model of polar liquids math- 
ematically. In following Section we discuss the bulk 
properties and the correlation functions of the liquid 
within the model. The results of the discussion can be 
used to calculate various asymptotic forms of interactions 
between point-like and macroscopic bodies of different 
shapes and charged states in a solute. At last, in Section 
4 we address the properties of hydration layers of molecu- 
lar dipoles on macroscopic surfaces with the emphasis on 
topological properties of hydrogen bonds networks and 
possible phase transitions. 



II. THE MODEL. 

Polar liquids are similar to ferroelectrics. Having this 
analogy in mind it is possible to develop a vector field 
theory in which the liquid is described by a local mean 
value of the molecular polarization vector 

s(r) = (d)/d (1) 

Here do is the static dipole momentum of a single 
molecule, and the averaging of molecular dipole momen- 
tums d is performed over a small but macroscopic volume 
of the liquid containing macroscopic number of molecules. 
Mathematical structure of theory of ferroelectrics 
is analogous to the phenomenological Landau-Ginzburg 
theory for ferromagnetics, which deals with the average 
value of atomic magnetic moment. All essential prop- 
erties of a polar liquid are the result of competition of 
two opposite effects, two types of interaction: the long- 
range electrostatics (dipole-dipole interaction) Q dd and 
the short-range intermolecular potential, £Ih, responsible 
for hydrogen bonds (H-bonds) formation. At room tem- 
peratures quantum-mechanical effects dominate in fin ■ 
As for fldd , it has pure classical nature. As the Earn- 
shaw's theorem says, the equilibrium state is impossible 
for the system of classical charges 0, 03| ■ It means that 
Qdd leads to chaotization of molecular dipole momenta 
orientations, and hence the ground state of a polar liquid 
is a disordered state at least on a large scale. At smaller 
distances however, the molecules tend to form maximal 
possible number of H-bonds, and thus the H-bond inter- 
actions order the molecular dipole moments. 

Below we employ the following expression for the free 
energy of the liquid: 

0(s(r)) = n H +n dd - J dVP Q sE ei (2) 

where E e is the external electric field, Po = nd n , and n 
is the density of molecules in a fluid. The general form 
of two first terms in l|2j is: 

{l H + n dd = ^Y, I dVdV's a {v)sp{v')M a p(r-r'). (3) 

a, (3 J 



Here a, = x,y, z enumerate the Cartesian components 
of the vector s. In a free liquid, E e = 0, Eqs.Q and © 
give the following expression for the correlation function: 

Q a p(v - r') = Mr) S/3 (r')) = TM~^(r - r'), (4) 

where T is the absolute temperature (in energy units) 
and the kernel is defined as: 

I dVxM-^v - rOAf 7/J (ri - r') = 5 a0 6{r - r'). (5) 

7 J 

The "matrix" M a/ 3 can be deduced from the correlation 
function taken, e.g. from molecular dynamics calcula- 
tions I n an isotropic liquid the correlation function 
has the form: 

Qa/3(r) = a(r)S a p + b(r)(f a fp - -S a fi), (6) 

where f = r/r. The functions a(r) and b(r) are related to 
the following correlation functions: a(r) = (s(0)s(r))/3, 
b(r) = 3(D(r) - o(r))/2, and D(r) = ((fs(0))(rs(r))). 
The functions a(r)and D(r) can be expressed with the 
help of the full pair correlation function g(r, 9\, 62, $) in- 
troduced in [13. Here 61, 82, </>iand 02 = <fri — <t>%) 
are the polar and azimuthal angles between vectors of 
molecular polarities, si = s(0), S2 = s(r), and the in- 
termolecular axis r. The correlation function g can be 
expanded in spherical harmonics: 

ff(r,0i,02,$) = 47r Y 9iihm(r)Y hm (6 U (t> 1 )Y l2 _ m (e2,<f>2)- 

(7) 

So that, gooo i r ) is a usual pair correlation function, which 
does not include information about molecular polar vec- 
tors. The correlation function Q gives the probability 
dW to find two molecules at a distance r with polariza- 
tion vectors close to si, S2: 

3ooo (r) 4tt 4tt ' 
Applying all the definitions above we find: 

a ( r ) = „ [2ffm(r) +3no(r)] /soooO): 



H r ) = 3 [ffuo(r) - ffin (r)] /goooM- 

The functions a and 6are unique characteristics of inter- 
molecular interactions. 

At sufficiently large scales the free energy functional 
has a few universal features. Since the dipole mo- 
ments of water molecules are ordered at small scales by 
hydrogen bonding, the mean molecular polarization vec- 
tor s(r) is a smooth continuous function. Therefore, the 
free energy of the liquid should allow expansion in powers 
of the molecular polarization gradients. That is why we 
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can use the Oseen's like general expression (see e.g. 
for the hydrogen bonding energy: 



ds 



2 dxp dx c 



V(s 2 



(8) 

where C, C are constants, and the scalar function V(s 2 ) 
takes into account short-distance part (other than dipole- 
dipole interactions) of the intermolecular interaction po- 
tential. The long-range part (dipole-dipole ) of inter- 
molecular interaction energy takes form: 



<>,;,; = / dV— E : 



p- 



(9) 



so that the full electric field in the liquid, E is the sum 
of the polarization and the external fields: E = Ep + E e . 
The polarization electric field Ep, satisfies the Poisson 
equation, having as a right-hand-side the density of po- 
larization charges pp = — divP, where P = PqS is the 
polarization of the liquid. 

In the bulk of a liquid in weak fields (E <C Pq) the vec- 
tor model defined by Eqs.{2J,(H| an d © can be studied in 
linearized approximation. First we expand the function 
V(s 2 ) in powers of s 2 : 



V(s 2 ) = ^ 2 (r), 



(10) 



where A is a dimensionless constant related to long-range 
correlation properties of the liquid (see below). Fourier 
transforming the model we find that: 



p2 

2 ^ 



s a Sf3 ({Ck 2 + A)5 a/3 + C'k a k ) + 



+iPo 22 k a a a <f> + ^2 



2A2 



k.a 



8tt 



where k a is the wavevector, 4> is the electrostatic poten- 
tial. The potential Sl ext describes the interaction of the 
linearized liquid with external objects. 



III. LONG RANGE INTERACTIONS 
BETWEEN SOLUTES. 

As seen from EqS.JHJ,® an d G3 polar liquid is nat- 
ur ally c haracterized by the two important scales: Lp = 
y/C/A and i?p = y / C/(A + Air). Indeed, on a molecular 
scale in small regions the dipole moments of the molecules 
are correlated and there exist pretty large electric fields: 
Ep ~ P ~ 1 10 7 V/cm. The size of such polarized region 
(a domain) is determined by the competition of dipole- 
dipole electrostatic and short-range hydrogen bonding 
forces: inside the domain fin ~ ^dd- Since fip ~ CPq R 
and Qdd ~ PqRq) where Rq is the size a polarized region, 



the scale Rd ~ \fC — 7A (the estimate C ~ 5 10~ 15 cto 2 
can be obtained from the value of the surface tension) 
can be considered as a typical size of domain (the ex- 
istence of such domains in a bulk of polar liquids was 
pointed in 20]). The second scale, Lp, includes ~ 100 
molecules and describes the correlations between such 
macroscopic domains and can be called as superdomain. 
It is the longest scale of intermolecular correlations in wa- 
ter (detailed consideration will be published elsewhere): 
at larger distances thermal fluctuations dominate. Super- 
domains determine a reaction of polar liquid on a weak 
uniform static electric field, i.e. its static dielectric con- 
stant. Experimental data for frequency dependence e(u>) 
are discussed in [SI US E3] ■ At room temperature the 
transition from the static value ~ 80 to the universal low 
value ~ 5 characterizing the internal rotational and elec- 
tronic degrees of freedom of individual water molecules 
is spread over a broad frequency range: 



10 5 < to < 



10V 



(11) 



Such sophisticated behaviour can be described with a 
help of a simple Debye model of polarization relax- 
ation (see pjj |2J], for example). In a weak electric 
field E the time evolution of liquid polarization fol- 
lows the linearized phenomenological kinetic equation 
P = — (P — (e — l)E/47r) /tq, where tq is the Debye re- 
laxation time. The solution gives the dielectric function 
e(u>) = (e — iiOTo) / (1 — iujtq) characterized indeed by a 
broad transient region with To ~ 10 _6 s. The latter quan- 
tity is much longer than the average life-time of a single 
hydrogen bond, ~ 10 _9 s, and thus should be related to 
water domains rearrangement processes. 

The model and 11011 is both non-local and non- 

linear. Full analysis of the solutions requires numerical 
simulations and is beyond the scope of this Letter. In 
what follows we confine ourselves to the solutions of lin- 
earized form of the equations. In spite of being quite a 
drastic simplification, this approach reveals the asymp- 
totic long-range interactions between solute objects and 
shows the importance of topological excitations, both in 
the bulk and on a surface of the liquid. 



A. Interaction between charged objects. Solvation 
energy of a point charge. 

Consider first the simplest case of interaction with a 
number of external point electric charges q a placed in 
positions r Q : 



n ext = J d 3 rp(r)(/)(r), 

where p(r) = J2a1 a ^( r ~ ra )- Then the free energy of 
the liquid is (setting C = for clarity): 



n = £>(k) 



, 47T 

e(fc)fc 2 



(12) 
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where 



e(k) 



R 2 D k 2 



Liipk 2 + 1 ' 



(13) 



is the effective dielectric constant at wavevector k, and 
e = L\jR 2 D = 1 + 4tt/A. The relation between the ma- 
terial constant A coming from the "equation of state" 
II10II and the macroscopic dielectric susceptibility e can 
be seen from the following observation: For a pair of 
point charges separated by a large distance R Eq. flSt 
gives the following expression for the interaction energy: 



Q(R) 



<M2 

R 



1 r. In ( R 

- + (l--)exp - — 
e e V R D 



(14) 



That is why at large distances R ^ Rd two charges 
interact as in a dielectric media characterized by dielec- 
tric constant e. In a polar liquid e> 1 and the constant 
A< 1 (for example at room temperature in water e w 80 
and A = 47r/(e- 1) « 0.16). 

The relation between A and e can also be established 
from a well known relation between the asymptotic be- 
havior of linear response functions of a liquid (correla- 
tors) and the dielectric constant (see e.g. [17ll25i p, Since 
A« 1, Rd <C Lt ~ 15A. Subtracting the vacuum en- 
ergy of a point charge from En. ltl2ll find the solvation free 
energy of a point charged object: 



n. 



k 



|p(k)| 



,47T 



1 



1 



1 



1 



2i? D 
(15) 

The result looks pretty much the same as classical Born 
solvation theory result. In fact, the solvation energy is 
localized in the electric domain of size k" 1 ~ Rd around 
the charge. In practical case Rd does not exceed much 
size of an ion (or a charged molecule). That is why the 
exact value of ionic Born radius Rb = ~q 2 /2Q S oiv ~ Rd 
depends on precise details of microscopic interactions of 
the ion and the surrounding water molecules. The results 
(1 14:1 and 11 51 come from the linearized model and hence 
are only valid for point-size charges q < qc = CPq 26] . 



B. Interactions of neutral objects. Typical types of 
liquid polarization. 

The interaction of macroscopic objects is not confined 
by electrostatic interaction of the charges. As demon- 
strated by direct numerical simulations, macroscopic ob- 
jects polarize the liquid 0,0,0. For example, next to a 
macroscopic surface the water molecules arrange them- 
selves such as their dipole moments d are parallel to the 
boundary (s 7^ 0, s _L n) with water molecules planes 
being perpendicular to the outer normal n to the surface 
of a body. The explanation is as follows. The bulk wa- 
ter molecule has N H w 4 hydrogen bonds with others. 
If the molecular polarization is orthogonal to the solute 
body, d || n, then Nh ~ 2 , and if d _L n, then Nh ~ 3. 




Figure 1: Polarization s of the liquid around a nearly spherical 
object. Forceless (a) and longitudinally polarized (b) liquid 
configurations. 



Therefore the liquid can lower its free energy by arrang- 
ing the molecular dipole moments along the hydrophobic 
surfaces (if no energy benefit can be extracted by hydro- 
gen bond formation with the solute atoms). This simple 
picture gives the following estimation for the surface ten- 
sion coefficient: a ~ (Eh /2)n~ 2 / 3 ~ PqRd, where Eh 
is the binding energy for one hydrogen bond. For water 
Eh ~ lOkJ/mol and a ~ 100 CGS, which is consistent 
with the macroscopic value a 70 CGS. 

In the presence of impurities the induced polarization 
of molecules propagates to a certain distance inside the 
bulk of the liquid and can cause a considerable interac- 
tion of hydrophobic objects at large distances. To analyze 
it we proceed to linear response study of polarization- 
polarization correlations in our model. At large distances 
the interaction of external objects does not depend on the 
precise structure of interacting bodies and takes univer- 
sal form. The asymptotic expression for the interaction 
potential can be obtained by observing the free energy 
change in a system of point (< Rd) chargeless impuri- 
ties located at the positions r a characterized by a single 
vector property j a . The simplest form of the interaction 
potential is 



^ext = / (j, s)dV, 



(16) 



where j(r) = J2 a iaS( r — r a ). For a specific shape of 
an interacting object the value and the direction of the 
vector j can not be found within the linear theory and 
requires microscopic derivation with the help of either 
complete model, or molecular dynamics calculation. 

In the linear response approximation the polarization 
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of the liquid is 



where 



G 



(8,8) 

a.8 



s a (r) = -J2 GapMMV exp(ikr), (17) 

k/3 



Po 1 f r 4irk a k (fc 2 C'/47r+l) 

On 



Ck 2 + A 



>a(3 



k 2 (C + C")k 2 + A + 4ir 
(18) 

The interaction of the two impurities separated by a large 
distance R = ri r 2 is given by 

n L (R)=- G ( :f{k) 3a] p = 

k,a,0 




b) 



E=0 




1 (ji,j2)-3(j 1 ,R)(j 2 ,R) 



p2 




A 2 eR 3 



This means that the two neutral impurities ill fit interact 
in unscreened fashion as two dipoles in a media with di- 
electric constant e. The dipole moment d arises from 
spontaneous polarization of water around the object: 
d = j/(APq). We call such liquid configurations as longi- 
tudinally polarized states, see e.g. FigQJ) for an example 
of water polarization force lines around a nearly spherical 
body. For such polarization type the non-vanishing po- 
larization charge exists in a liquid: pp = — Podivs ^ 0. 
The induced dipole moment around a microscopic ob- 
ject of the size Rq < Rd can be estimated as follows: 
at small distances r < Rd (k — > oo), the second term 
in the correlation function iflSll can be neglected and 
s ~ j/ {PqCt). The linear theory breaks on the surface of 
the object, r = Rq, where s ~ 1. Therefore j ~ P 2 CRq, 
d ~ PqCRqJA and the energy of a pair of chargeless ob- 
jects can be estimated as 



n L (R) 



C 2 R 2 a P 2 
AR 3 ' 



(19) 



This interaction is long-range and is completely due to 
the dipole-dipole interaction of the water molecules. 

Longitudinal polarization of the liquid also contributes 
to the solvation energy of a single impurity: 

ftsoiv = -J2 ^jpG { ^f (k) = n + n A , (20) 

k,a,/3 

where fio comes from the first term of Eq. I|18|) , formally 
diverges and can be estimated as fio ~ j 2 Pok max /C, 
where k max ~ R^ 1 is related to the size of the impurity. 
The contribution from the second term is finite and is 
in fact the energy ^ of the polarization electric field: 
flA ~ P 2 C 3 / 2 = 50 kJ/mol. Since ft ~ CRqP$, Vt A > 
Slo for impurities of sufficiently small size Rq < Rd, the 
quantity £Ia is large and can be called as the activation 
energy of the longitudinally polarized configuration. 

Large value of the activation energy is due to appear- 
ance of strong electric fields next to polarized bodies. In 



Figure 2: Polarization of a liquid near a cylindrical hydropho- 
bic surface. Forceless (a) and longitudinally polarized (b) liq- 
uid configurations. 



fact there is a wide class of liquid configurations with 
appreciably lower energies due to the absence of polar- 
ization charge. This can happen in a system of neutral 
bodies if a special, the "forceless" (FL), state of molecular 
dipoles is chosen: 



Vs 



0, E = 0. 



(21) 



Two examples of forceless water configurations around a 
sphere and a cylinder are shown on Figs. The po- 

larization vector in such states is similar to magnetic field 
in magnetostatics, therefore the asymptotic form of the 
interaction potential with with external point size water 
polarizing impurities can be specified using the following 
form of the interaction potential: 



^ext — 



y*(J,V x s)dV. 



(22) 



where J = ^ a <5(r — r a )J a , r Q is the position of an im- 
purity, and J a is a vector property of an impurity . The 
vector J a depends on the details of the surface-water in- 
teractions and requires full microscopic calculation for its 
determination. The energy of interacting bodies is then: 

= ^ ] ^ {fytapvtpii'v' k^J v k^i J v i , 

k,a,0... 

where e Q /3 7 is the totally antisymmetric tensor. Using 
Eq. I|18|) we find that the contribution of the second term 
vanishes. At intermediate distances R < Lp 

Qfl{K) = A^pJCR 3 (23) 

At larger distances R > Lt the interaction vanishes 
exponentially: 



o ro \ (Ji,R)(J 2 ,R) - (Ji, 3 2 ) R/L 



AtiP 2 CRL 2 t 



(24) 
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The estimation for the vector J proceeds as follows: at 
small distances s ~ J/(PqCR 2 ) and has to be s ~ 1 
at the surface R ~ i?o of the body. Hence, Ji,2 = 
AttPqCRI 2 n i,2, where 111,2 are the unit vectors con- 
nected with the vorticity of polarization (see Figs, and 
Et) by the "right-screw rule", R\ : 2 are the characteristic 
radii of the objects. The interaction of a pair of charge- 
less impurities in a forceless water configuration at small 
112311 and large Il24t distances reads, correspondingly: 

n™ ~ ^ (25) 

riT/A p2 

RL^~ ' ^ ' 

The exponential decay of the interaction at large dis- 
tances is not surprising since forceless configurations are 
also chargeless and the net dipole momentum of the sys- 
tem is identically zero. Note, that for a forceless configu- 
ration the polarization of the liquid is present at distances 
A ~ Lt from the solute surface. 

Eqs. p^f . lj2"fijl universally characterize forceless inter- 
actions of small neutral objects. In linear approximation 
l|17H the polarization vector around a pair of solutes lo- 
cated at points pi and p2 is given by 

s(r) = s FL = Rl{n 1 y.r l )f{r l ) + R 2 2 {n 2 y.r 2 )f{r 2 ), (27) 

where f(r) = r~ 2 (l + r/Lx) exp(— r/Lx). Here ri,2 — 
r — pi, 2, ri,2 = 1*1,2/7*1,2, and R\ : 2 are the characteristic 
radial dimensions of the solute particles. Eq. il2?jl is valid 
far from the objects (7*1,2 -^1,2)- If the solute bodies 
carry charges, than additional induced polarization con- 
tribution appears: 

s =s FL + sp, (28) 



sp ~ -8ig{n/R D )ri - 8 2 g{r 2 / R D )v 2 , (29) 

where g(x) = x~ 2 [1 - (1 + x) exp (—x)], 5i t 2 = qx^lqc- 
The water polarization is qualitatively represented on the 
Fig. 03 Water polarization of this type was obtained in 
MD simulations ^3] and called as "the dipole-bridge be- 
tween biomolecules". The vortex-like structures of polar- 
ization were observed that corresponds to the term sp^ 
in Jl§. 

At intermediate distances between the interacting bod- 
ies both forceless (FL) and longitudinally polarized (L) 
liquid configurations have very similar interaction prop- 
erties. Both and Qfl are dipole-dipole interactions 
of polarized molecular dipoles. For sufficiently small im- 
purities we have Q.fl/^l ~ Rq/Rd ^ 1- This means 
that the interactions of small longitudinally polarized 
impurities is stronger. On the other hand, the second 
term of Eq. ltl8)l does not also contribute to the solva- 
tion energy of a single forceless impurity, which means 




Figure 3: Polarization of a liquid near two hydrophobic 
charged solutes. 

that the forceless configurations does not have the large 
activation energy contribution Qa S> T in its free en- 
ergy and hence, at least at large distances, thermody- 
namically is more probable. At intermediate distances, 
R £ L T {Rl/R D L T ) 1/3 < L T , fl L > Q A and sponta- 
neously polarized liquid configuration with anti-parallel 
j possesses lower energy and may become favorable. This 
means that anisotropic neutral molecules in a polar liquid 
can acquire fairly strong orientation dependent interac- 
tion caused by induced dipole moment of the surrounding 
water molecules. 



C. Interaction of macroscopic objects: parallel 
hydrophobic cylinders. 

The interaction potentials l|19Jl and l!24t give only the 
energies of point-like objects of size Rq <C Rd- In re- 
alistic calculations Rq > Rd (small organic molecules) 
and may even be larger (proteins, lipids etc.). For larger 
objects, a more detailed calculation has to be done. 
Consider first two parallel cylinders of the same radii 
i?i,2 = Rq, heights H\$ — H ^> Rq, placed at a distance 
L <C H from one another. The separation between the 
objects is large, L ^> Rq, but the size Rq is now arbitrar- 
ily related to the water domain size Rd- At distances 
large enough from the cylinder cores the superposition 
principle holds: 

s(r) wsi(r)+sj(r), (30) 

where Sj, are cylindrically symmetric forceless solutions 
of model around isolated cylinders (see FigEk): 

Bi(r) = Xi—(3i$i. (31) 
n 

Here Ai,2 = ±1 stands for two possible directions of wa- 
ter dipoles around a cylinder (the "topological charges"), 
Ti are the distances from the point r to the axis of 
the correspondent cylinder. The coefficients: (3\.2 ~ 
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1/ [1 + jRi^/Rd], where 7 ~ 1 depends on the be- 
haviour of the function V(s 2 ) for s ~ 1. Formula 1(31")) 
holds for Ti <C Lt- Since Vs « 0, the polarization elec- 
tric field vanishes. Substituting directly the solution l-'iOli 
into the energy density <|2J we find the following expres- 
sion for the interaction potential: 



n FL (L) « 2ttAP 2 \ 1 \ 2 (3 1 i3 2 R 2 HK 



L 



(32) 



where K is the Macdonald function. For small cylin- 
ders, R <C Rd, (3 « 1 and the result ll32)l coincides with 
that one can find by integrating the potential 12411 over 
the cylinders in the linear model. At sufficiently small 
distance between the cylinders the interaction exceeds 
the activation energy fi^ and the longitudinally polar- 
ized configuration settles (see Fig. Eb)- Here Vs =^ 
close to the ends of the cylinders, so that the energy of 
each of the cylinders contains ~ 2ttRI^/CttPq . The 
superposition principle still holds and the energy of the 
cylinders reads: 



Cl L (L) » 2ttCP 2 ^P 1 (3 2 HK 
A1A2 



L 

Li ' 



(33) 



Here Ai,2 = log(Lr/-Ri,2)- The result is again similar to 
that one can find by integrating the interaction poten- 
tial for a spontaneously polarized configuration l|19JI over 
entire cylinders length. 

Egs- EOjl and l(3T|l describe an approximate forceless so- 
lution. Let us estimate the leading correction fl' to the 
interaction energy due to the dipole-dipole interactions 
of polarization charges (divs ^ 0). Consider first two 
cylinders at a distance L < Lr- In this case the po- 
larization charges are concentrated around the cylinders 
surface, Qa.b — —Qc.d ~ PqRo in the sectors A — D 
of Figure El The dipole moment D ~ PqRqL per unit 
length along the cylinders gives the following estimation 
for the interaction energy 



Q' 



e{k ~ L- 1 ) 



dz\dz'_ 



[(DiD 2 )-3(nD!)(nD 2 )] 



1 12 



2H 



e{k~L- l )L- 



[(DiD 2 ) - 2D lx D 2x ] (34) 



The integration in l(3"4"jl is performed along the axes zi ; 2 
of the cylinders. The vector r 12 connects the points at 
the axes, ri2 = y ( z i — z z) 2 + -^ 2 ; an d 11 = r i2Ai2- For 
a pair of cylinders of the same size Di = T> 2 = D and 



n' 



P 2 HR 2 (L 2 + L 2 T ) 



e(L 



2 + R 2 D ) 



(35) 



At larger distances, L > Lt, the solution l(3"l")l acquires 
an extra exponential factor exp(— r/Lx) and 

p2 ZTE>2 

Sl'~ exp(-2£/L T ). 



Comparing this result with Eas, l|32JI and ll-'iot we find: 

n'(L)/il FL (L) ~ exp(-£/L T ) < 1. 

This means that for sufficiently large distances the po- 
larization charges provide negligible corrections to the 
energy of a forceless configuration. 



IV. SOLVENT-SOLUTE INTERFACES. PHASE 
TRANSITIONS ON BOUNDARIES. 

Consider now the structure of molecular polarization 
near a plane liquid-to- vacuum or liquid-to-a body surface 
interface. Such a boundary is hydrophobic and polarizes 
the liquid along the boundary plane T. The polarization 
extends to the bulk of the liquid on a distance scale A ~ 
Lt or A ~ Rp depending on electrostatic properties of 
the boundary (see below). Since at larger distances from 
the boundary the polarization of the liquid disappears, 
the polarization itself is an effectively two-dimensional 
vector field. The Hamiltonian for the surface polarization 
field s ]| can be obtained by integration of {01 inside the 
bulk of the liquid: 



1 



1 



M / df (yey + -K / df (V • S) 0(r) (36) 



0(r) 



df 



,(V'-S') 



(37) 



Here 0(x,y) is the angle characterizing the direction of 
the unit vector S = sy/sy = (cos 9, sin 6), r = (x,y) 
are the coordinates along the surface, V = (d/dx, d/dy), 
4>(r) is the electric potential at the point r. The constants 
M ~ CP 2 As jj, K ~ AM . Both terms in Eq. JUJ) 
originate from Ojy and Sldd of Eq. |2j , respectively. Note 
that 



V-S 



sinfi 1 ■ 9 X + cos 9 ■ 9 y , (38) 



where 9 X = 89 /dx, 9 y = 89 jdy. Minimization of the 
functional lj3fil) with respect to the variations of 9(x, y) 
gives the following "equilibrium condition": 



MA9(x,y) + K(-sm9 



cos f 



/) = (39) 



The simplest solutions of the selfconsistent equations 
(l39l) . (JHHI is the uniform polarization: 9(x,y) 

cons' and 



const 



(40) 



A more sophisticated solution describes a vortex state: 
9(x,y) = m arctan(y/x) of topological charge q = m = 
0, ±1, ±2, .... Both mentioned solutions are forceless con- 
figurations and therefore A ~ L T . The vortex core size 
is also ~ L T . Fig0 shows a more complicated polariza- 
tion configuration, which is is a vortex-antivortex pair. 
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Figure 4: Polarization configuration S(x,y) of surface water 
"vortex- ant i vortex pair". 




Figure 5: Excitation of water with one vortex- ant i vortex pair 
pinned on a surface of hydrophobic solute. 

In this case V-s ^ and thus A ~ Rd for all essential 
configurations with vortex-antivortex pairs. 

Let us neglect first the dipole-dipole interaction term 
Qdd in l(3fijl . In this case En. ll36|l coincides with the 
Hamiltonian of XY 2D model. As it is well known from 
the field theory |27l l2^|. thermal fluctuations of the po- 
larization field iKffiJl can be mapped on to the thermal 
dynamics of a gas of interacting vortices residing on the 
surface. The thermal state of the liquid mainly consists 
of vortex-antivortex pairs of minimum charge \q\ = 1. 
The energy of vortex-antivortex configuration at large 
distance 

L > A (41) 

between their cores is given by 

Q(L) w ttx +vrMlog(£/A) , (42) 

where fli ~ M, We observe that the "hydrogen-bonding" 
term Q.h leads to appearance of the logarithmic attrac- 
tion of vortex and antivortex. At temperature 

Tbkt = \kM ~ ttCP 2 A, (43) 



the topological Berezinskii-Kosterlitz-Thouless transition 
[22], [23] occurs: bound vortex-antivortex pairs dissociate 
and form the vortex "plasma". According to Onsager ar- 
guments [11113,113 e ~ p o r d/ t - On a hydrophobic 
surface s?, ~ 1, and 

Tbkt ~ eT » T. (44) 

This means that at any reasonable temperature the po- 
larization field remains ordered (the molecular dipoles are 
correlated at large distances, see below). 

The long-range dipole-dipole interaction described by 
the second term in Ea . essentially alters the character 
of BKT transition (221 . As is clear from FigQ]the polar- 
ization vector si produced by the left vortex encounters 
other vortex at the region next to the point C and thus 
induces the polarization density pp{C) > 0. At large 
distances fill) only one scale L is essential: \rc — i\d| ~ 
|rc — tb\ ~ L, therefore pp{C) ~ s\\Pq/L. The polar- 
ization charges in each of the sectors equal: Qc = Qd = 
—Qa = —Qb ~ pp{C)XL 2 ~ s\\Pq\L. Roughly speak- 
ing the polarization charge pattern can be approximated 
by the charge distributions of two equal downward ori- 
ented dipole moments D ~ QcL ~ s^PqXL 2 placed at a 
distance ~ L one from another. Simple estimations give 
tidd ~ D 2 /L 3 - KL and 

f2(£)psn 1 + 7 rMlog(L/A)+ 0KL, (45) 

with /3 ~ 1, instead of (@2j. In other words the dipole- 
dipole interaction of vortex with antivortex gives addi- 
tional attraction term with a linear dependence on L. 
The following remark is necessary to explain this inter- 
action. Indeed, the charge distribution on the left side 
of Fig. 01 is a "mirror reflection" of that on the right 
side. Such charge distributions repel in vacuum. In the 
case of a polar liquid interface the polarization charges 
Qc increase as L increases so that Qc + Qb = at any 
L. It means that in order to move apart a vortex and 
antivortex pair, one should produce additional work to 
polarize the water molecules. The Mayer-Schwabl phase 
transition occurs at 

T = T MS = 4T B kt- (46) 

The phase transition is analogous to that for quark-gluon 
plasma formation [13]: unbound vortices are formed at 
T > Tms- To establish this fact one has to go beyond the 
approximation of pair interaction of vortices jjH - Ac- 
cording to lH5)) the typical dimension of vortex-antivortex 
pair is: Lp ~ T/K. The surface concentration of pairs is: 
np ~ A~ 2 exp(— Qi/T). At npLp ~ 1 we have a plasma- 
like quasineutral system, in which topological charges of 
vortices and antivotices are mutually compensated. Us- 
ing our estimations for the parameters K, Oj and A we 
conclude that at T > Tms the interaction J45fl between 
two vortices is modified and the phase transition occurs. 
In fact, Eq. lllfi)) between describes the interaction of vor- 
tices only at T < Tms ■ At higher temperatures above 
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the phase transition point we can use a set of standard 
Debye arguments as in linear approximation of screening 
in plasma (see e.g. [19]). and we find the following mean- 
field effective interaction between vortices of Qt placed 
at a distance L from each other : 



fi(L) 



Qt 



dpp 



pJo(pL) 



(p 2 + K 2 )p + K 2 pi 



Here k = v / 2Trn P M/T ~ 1/A and pi = K/(wM) ~ 1/A. 
Note that k ~ p x at T ~ T MS . At T ^> T MS n ^> p 1: 
and at small distances, L <C 1/pi, the the interaction is 
logarithmic: fi(L) » ttMKq(kL) + (3K L , which recovers 
Eq. l(4*5l) at sufficiently small distances L -C 1/k. At 
large distances L ^> 1/pi the screening switches on and 
the interaction decays as Q,(L) w Qt/(k 2 PiL 3 ) ~ 1/i 3 . 

Eqs. l|l4Tl and ijlij) show clearly that the phase tran- 
sition can not be observed on a hydrophobic surface. 
Nevertheless it may occur on hydrophilic surfaces, where 
the energy of the liquid can be further lowered by extra- 
hydrogen bonds involving the atoms on the surface of the 
solute. In a place of such polar contact the vector s is al- 
most perpendicular to the liquid boundary: s| -C 1. The 
effective Hamiltonian for the molecular dipole field on a 
hydrophilic surface still has the form Il3ti|) . though the 
phase transition temperature is lower: Tms ~ Ci^o^f j 
Tms/T ~ esjj. The appearance of the additional small 
factor s| can decrease Tms drastically: the transition 
temperature may be shifted to the room temperatures 
range already on a hydrophilic surface with sii ~ 0.1. 

At T > Tms the vortices dissociate and form vortex- 
antivortex plasma, so that the order parameter s fluc- 
tuates strongly. One may say that at small tempera- 
tures (or on a highly hydrophobic boundary) the molec- 
ular dipoles fluctuations are small and the molecular di- 
rections field is highly ordered on large scales. In other 
words global network of hydrogen bonds exist on the sur- 
face of the solute body. As the polar interactions of the 
solute object and the surrounding liquid increase, the the 
phase transition temperature gets lower and the global 
polarization order, and hence the hydrogen bonds net- 
work, disintegrates. Moreover the hydrophilic regions on 
the surface of the solute may pin the vortices, giving rise 
to a peculiar sophisticated polarization shapes (see Fig. 
© and the discussion in |Kl|V 

We argue that the above discussed phase transition 
is directly related to the disruption of the global hydro- 
gen bonds network in the course of 2D percolation phase 
transition observed in the Molecular Dynamics simula- 
tion of hydration water absorbed on the surface of a so- 
lute [ill fia. flT^ I . In our model this result can be explained 
with the help of classical field theory considerations. At 
T > Tms a solute boundary hosts a gas of independent 
vortices of different "circulation" resulting in fast decay 
of correlation between dipole moments of distant water 
molecules: D(r) = (d(O)d(r)) - exp(-r/r c ) [j^. At 
T < T M s, the vortices annihilate the and long-range cor- 
relations set in: the ferromagnetic-like ordering of hy- 




Figure 6: Ferroelectric wave propagates in a film of hydration 
water. 



dration water molecules dipole moments takes place. It 
means that in the considered case of polar liquids the 
boundary layer of liquid is a ferroelectric film. Such cor- 
related thermal states can be seen as global network of 
hydrogen bonds at the surface. Note that in our model 
the value of transition temperature Tms is proportional 
to the width of polarized layer A, which is in accordance 
with the dependence of perc olation transition tempera- 
tures obtained in H3.ll^ as a function the hydration 
shell width D (of course, it is reasonable to suppose that 
D ~ A). 

The appearance of ferroelectric order can be used as a 
justification for the mean field picture employed through- 
out the Letter. This fact does not contradict with the 
Peierls-Mermin theorem [34, 35] on the absence of order- 
ing in a 2D systems due to new effect: the long-range 
dipole-dipole term in l.ifill . To illustrate this point con- 
sider collective vibrations of molecular dipole moments in 
a ferroelectric hydration water layer (FigEJ with molec- 
ular dipoles arranged along x axis. The kinetic energy 
of a water molecule is I6 2 /2, where I is the molecular 
moment of inertia. The dispersion law oj = to(p) of the 
ferroelectric wave of a wave vector p can be established 
by the minimization of the action S — J Ldt, where the 
Lagrangian L reads 



L 



df l -N9 2 - fi s . 



Here N = An J and the interaction potential fls is defined 
by l(3fijl . The equation of motion 

-N6 + MA9(x, y) + K(- sin 9 ■ <j> x + cos 6 ■ <f> y ) = 

can be expanded in powers of small \9\ <§; 1 so that the 
dispersion relation is: 



IMp 2 2irK . 2 

w(P) = \ 1 psin a. 

yv ' N N 



(47) 



Here a is the angle between vector p and x axis. Consider 
the Peierls-Mermin argument applied for a 2D atomic 
crystalline lattice. Mean squared displacements of the 
atoms from the equilibrium positions in the lattice are 
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proportional to the integral J = J d 2 p/u> 2 (p). Short- 
range interactions of atoms give the phonon excitations 
with uj ~ p. In this case J diverges at small values 
of p, which corresponds to fluctuations with large wave 
lengths. The integral with the dispersion relation I|47J1 
does not have the singularity at small p: 

t*2ir />oo />oo 

J~ / da dp/(Mp + 2nK sin 2 a) ~ / dp/^/p, 
Jo Jo Jo 

(48) 

which means that the long-range dipole-dipole together 
with the short-range hydrogen bonding interaction both 
stabilize the ferroelectric ordering of the film in agree- 
ment with Interestingly both Qh and Qdd are 
equally important in this phenomenon: long range or- 
der does not exist both in the film with the short and the 
long range interactions separately: the integral diverges 
either if M — (no short range hydrogen bonds in the 
model) and if if = (no dipole-dipole interaction). This 
somewhat surprising observation means that the short- 
range behavior of intermolecular interactions determine 
the character of order at large scales. Physically it is a 
consequence of the Earnshaw's instability of the system 
of classical dipoles, which is stabilized by the H-bonds. 
Similar findings are reported in [sfj l3?l Isfl. f3Si | where the 
2D system of classical dipoles were considered on a few 
different types of lattices. It was shown that the ground 
state of the lattice system can switch from a ferromag- 
netic to antiferromagnetic state depending on the lattice 
shape. 

V. CONCLUDING REMARKS. 

Finally we compare the results of our vector model 
to the density functional theory, which is another well 
known approach to hydrophobic interactions calcula- 
tions. It was pioneered in earlier works of van Der 
Waals [40], works well for simple, non polar liquids with 
e ~ 1, does not focus on electric properties and uses the 
free energy functional expressed in terms of the mean 
n(r) = (p(r)), the total p(r) density of the liquid, and 
the density fluctuations w(r) ((w(r)) = 0) 0, 0,0,111 E3- 
The free energy functional for the mean density is a non- 
linear function of n(r) and is characterized by the healing 
length a mac ~ 7A for water. The fluctuations w(r) are 
described by the pair correlation function taken from ex- 
periments or molecular dynamics simulations and is char- 
acterized by the second, microscopic scale, a m ic ~ 1A. If 
the interaction between the molecules of a liquid is purely 
short range, the correlation function for the fluctuations 
decays exponentially at large distances: 

n I tn\ i \\ cx P(- r /Qw) 
G = (w(O)w(r)) . 

The interaction potential fl of a pair of point-like ob- 
jects is proportional to the correlation function, O(r) ~ 
G(r) ~ exp(— r/a m i C ) and, due to the scalar nature of 



underlying theory, is not orientation dependent. Techni- 
cally speaking the model suggested in this paper also de- 
scribes a model system with two characteristic scales, say, 
a-mac ~ Lt and a m i C ~ Rjj and thus can share many pre- 
dictions with the density functional models. On the other 
hand, as was shown above the long-range dipole-dipole 
interaction between the molecules of the liquid leads to 
a long-range power-law tail in the correlation function 

{s a sp) ~ -3 , r -> 00. 

This property is inherently related with the large dielec- 
tric constant and is directly responsible for the strong 
long-range orientation-dependent interactions described 
in this Article. Strictly speaking the liquid in our model 
is incompressible, the assumption obviously fails at r ~ 
a m i C . The complete model should include both fluctua- 
tions of the liquid density (scalar) and the polarization 
(vector) . One example of such a model of a polar liquid 
was studied in |2q . Molecules were considered as freely 
oriented, like in the Onzager's theory 0, H3, OH • The 
approach proved to successfully reproduce the results of 
molecular dynamics simulations for solvated ions. An- 
other route to include the effects of liquid polarization 
and electrostatic forces in the scalar model is outlined 
in 0| . The advantage of vector models is clearly seen 
when the effects of the molecular polarization need to be 
taken into account, e.g. in a case of charged solutes or 
for liquid boundaries and interfaces. No scalar model can 
be directly applied to hydrogen bonds network formation 
studies and topological excitations and interactions both 
in the bulk and on the solute boundaries. 

The proposed vector model is both physically rich, 
relatively simple, and is not much more computation- 
ally demanding than continuous solvation models used 
in computational biophysics (see, e.g. 0, 0|). In its 
linearized form it reproduces a few very important fea- 
tures of molecular dynamics simulations, such as sponta- 
neous liquid polarization next to macroscopic surfaces, 
the vortex-like structures in the molecular dipoles ar- 
rangement near soluted bodies. Strong interaction be- 
tween the liquid molecules selects forceless liquids config- 
urations with vanishing polarization electric field. This 
leads to a strong, long-range, orientation dependent in- 
teraction of solvated bodies. The model may serve well 
as a continuous water model to describe different oper- 
ational states of nanodevices in aqueous environments, 
as well as membrane charge transfer and ion channels, 
where both the hydrophobic and charged boundaries are 
equally important. 

Fast solvation models are crucial in drug discovery and 
biomolecules modelling. The quality of solvation models 
is a limiting factor determining the accuracy of the calcu- 
lations in molecular docking and virtual screening appli- 
cations. The proposed model is employed in QUANTUM 
drug discovery software for very accurate binding affin- 
ity predictions for protein- ligand complexes QUAN- 
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TUM uses the continuous vector water model (2J for sol- states of biomacromolecules, pKa calculations, protein 
vation energy predictions, determination of protonation structure refinement (including mutagenesis) etc. 
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